MODULE constants
  IMPLICIT NONE
  REAL*8,PARAMETER:: c_pi=3.14159265358979323846d0
  !REAL*8,PARAMETER:: c_pi=4.d0*datan(1.d0)
  REAL*8,PARAMETER:: c_try=1.23456789123456789123456d0
  REAL*8,PARAMETER:: c_na = 6.02214d+23 ! Avogadro's number
  REAL*8,PARAMETER:: c_h = 6.62606896d-34 ! Planck's constant (J*s)
  REAL*8,PARAMETER:: c_c = 2.99792458d+10 ! Speed of light in vacuum (cm/s)
  REAL*8,PARAMETER:: c_h2wn = 2.194746d+05 ! Hartree to cm^-1
  REAL*8,PARAMETER:: c_amu = 1.660538782d-27 ! = 1/_na/1000
  REAL*8,PARAMETER:: c_gamma = 4.d0 * c_pi * c_pi * c_c * 1.0d-23 / c_h / c_na;
  REAL*8,PARAMETER:: c_bohr=5.291772108d-11!bohr radius in meters =
  REAL*8,PARAMETER:: c_angs2bohr=1d-10/c_bohr
  REAL*8,PARAMETER:: c_emass=9.10938215d-31!electron mass in kg
  REAL*8,PARAMETER:: c_amu2emass=c_amu/c_emass
END MODULE constants

MODULE global
  IMPLICIT NONE
  INTEGER maxcell,maxk,Ndof,Ndim,Natom,Ncell,N_k,nMR
  INTEGER, ALLOCATABLE::maxbasis(:)
  REAL*8, ALLOCATABLE ::mass(:),k_freq(:,:),hrmfreq(:),NC(:,:),coef(:,:)
  COMPLEX*16, ALLOCATABLE ::k_coef(:,:,:),cck_coef(:,:,:)
  REAL*8  :: geo(100)
END MODULE global

MODULE forceconstants
  IMPLICIT NONE

  REAL*8,ALLOCATABLE                    :: x_grad(:)
  COMPLEX*16,ALLOCATABLE              :: q_grad(:,:)
  REAL*8,ALLOCATABLE                :: x_harm(:,:,:)
  COMPLEX*16,ALLOCATABLE          :: q_harm(:,:,:,:)
  !  real*8,allocatable            :: x_cube(:,:,:,:,:)
  !  complex*16,allocatable     :: q_cube(:,:,:,:,:,:)
  !  real*8,allocatable       :: x_quar(:,:,:,:,:,:,:)
  !  complex*16,allocatable :: q_quar(:,:,:,:,:,:,:,:)

  !variables for nMR QFF, n=1,2,3,4
  CHARACTER*100 :: title1MR,title2MR,title3MR
  REAL*8:: E0
  REAL*8,ALLOCATABLE :: Gi(:),Hii(:),Ciii(:),Qiiii(:)
  REAL*8,ALLOCATABLE :: Hij(:,:),Ciij(:,:),Qiijj(:,:),Qiiij(:,:)
  REAL*8,ALLOCATABLE :: Cijk(:,:,:),Qiijk(:,:,:)
  REAL*8,ALLOCATABLE :: Qijkm(:,:,:,:)
  ! real*8 :: Gi(Ndof),Hii(Ndof),Ciii(Ndof),Qiiii(Ndof)
  ! real*8 :: Hij(Ndof,Ndof),Ciij(Ndof,Ndof),Qiijj(Ndof,Ndof),Qiiij(Ndof,Ndof)
  ! real*8 :: Cijk(Ndof,Ndof,Ndof),Qiijk(Ndof,Ndof,Ndof)
  ! real*8 :: Qijkm(Ndof,Ndof,Ndof,Ndof)

END MODULE forceconstants

MODULE integral
  IMPLICIT NONE
CONTAINS

  REAL*8 FUNCTION q_integral(iexc,m,power,diff)
    ! <n|Q_m^power|n-diff>
    USE constants,ONLY:c_gamma
    USE Global
    IMPLICIT NONE
    INTEGER ::iexc,m,power,diff
    REAL*8::n,n_2,n_3
    REAL*8::gamma,gamma_2,gamma_3
    ! <n|Q_m^power|n-diff>
    n=dfloat(iexc-1)
    n_2 = n * n;
    n_3 = n * n * n;
    gamma = c_gamma * hrmfreq(m)
    gamma_2 = gamma * gamma
    gamma_3 = gamma * gamma * gamma
    q_integral=0.d0

    SELECT CASE (power)

    CASE (1)
       IF (diff == 1) THEN
          q_integral=DSQRT(n * 0.5d0 / gamma);
       END IF
    CASE (2)
       IF (diff == 0) THEN
          q_integral=(n + 0.5d0) / gamma;
       ELSE IF (diff == 2)THEN
          q_integral=0.5d0 * DSQRT(n * (n - 1.d0)) / gamma;
       ENDIF
    CASE (3)
       IF (diff == 1)THEN
          q_integral=3.d0 * DSQRT(n_3 / (8.d0 * gamma_3));
       ELSE IF (diff == 3) THEN
          q_integral=DSQRT(0.125d0 * n * (n - 1.d0) * (n - 2.d0) / gamma_3);
       ENDIF
    CASE (4)
       IF (diff == 0)THEN
          q_integral=0.75d0 * (n_2 * 2.d0 + n * 2.d0 + 1.d0) / gamma_2;
       ELSE IF (diff == 2)THEN
          q_integral=0.5d0 * (2.d0 * n - 1.d0) * DSQRT(n * (n - 1.d0)) / gamma_2;
       ELSE IF (diff == 4)THEN
          q_integral=0.25d0 * dsqrt(n * (n - 1.d0) * (n - 2.d0) * (n - 3.d0)) / gamma_2;
       ENDIF
    END SELECT
  END FUNCTION q_integral

  REAL*8 FUNCTION kin_integral(level,m,diff)
    USE constants,ONLY:c_gamma
    USE Global
    IMPLICIT NONE
    INTEGER ::iexc,m,power,diff,level,n,n_2,n_3
    REAL*8::gamma,gamma_2,gamma_3
    ! <n|del^2/del(q)^2|n-diff>
    n=level
    n_2 = n * n;
    n_3 = n * n * n;
    gamma = c_gamma * hrmfreq(m)
    gamma_2 = gamma * gamma
    gamma_3 = gamma * gamma * gamma
    IF (diff == 0)THEN
       kin_integral=0.5 * gamma * (n + 0.5)
    ELSE IF (diff == 2)THEN
       kin_integral=0.25 * gamma * SQRT((n + 1.0)*dfloat(n) )
    ELSE
       kin_integral=0.0
    ENDIF
  END FUNCTION kin_integral
END MODULE integral

MODULE integral2
  USE integral
  IMPLICIT NONE
CONTAINS

  REAL*8 FUNCTION fullIntegral(m,power)
    USE global
    IMPLICIT NONE
    INTEGER ::m,power,n1,n2
    fullIntegral = 0.d0;
    DO n1 = 1,maxbasis(m)
       DO n2 = 1,maxbasis(m)
          fullintegral = fullintegral + coef(m,n1) * coef(m,n2) * q_integral(MAX(n1, n2), m, power, ABS(n2-n1))
       ENDDO
    ENDDO
  END FUNCTION fullintegral

END MODULE integral2

MODULE qVSCFintegrals
  USE global,ONLY:ndof
  USE forceconstants
  USE integral
  IMPLICIT NONE
CONTAINS
  REAL*8 FUNCTION uzero(mode)
    IMPLICIT NONE
    INTEGER::m1,m2,mode
    uzero=0.0
    DO m1=1,ndof
       IF (m1 .EQ. mode)CYCLE
       print*,'m1= ',m1
       uzero= uzero + hii(m1) * q_integral(1, m1, 2, 0)
       print*,'uzero= ',uzero
       DO m2=1,ndof
          IF (m2 .EQ. mode .OR. m1 .EQ. m2)CYCLE
          uzero= uzero + 0.125d0 * qiijj(m1,m2) * q_integral(1, m1, 2, 0)* q_integral(1, m2, 2, 0)
       ENDDO
    ENDDO
  END FUNCTION uzero
  REAL*8 FUNCTION utwo(mode)
    USE constants
    IMPLICIT NONE
    INTEGER::m,mode
    utwo=0.0
    DO m=1,ndof
       IF (m .EQ. mode)CYCLE
       utwo= utwo + 0.5d0 *qiijj(mode,m) * q_integral(1, m, 2, 0)
    ENDDO
    utwo= utwo + hii(mode)/c_angs2bohr/c_angs2bohr/c_amu2emass*c_h2wn*2.d0
  END FUNCTION utwo

  REAL*8 FUNCTION kinetic(level,mode)
    USE constants,ONLY:c_gamma
    USE Global
    IMPLICIT NONE
    INTEGER ::iexc,m,power,diff,level,n,mode
    REAL*8::gamma
    ! <n|del^2/del(q)^2|n-diff>
    n=level-1
    kinetic=0.d0
    DO m=1,ndof
       IF (m .EQ. mode)CYCLE
       gamma = c_gamma * hrmfreq(m)
       kinetic= kinetic+0.5 * gamma * (n + 0.5)
    ENDDO
    RETURN
  END FUNCTION kinetic
END MODULE qVSCFintegrals

MODULE yazar
  IMPLICIT NONE
CONTAINS
  SUBROUTINE writematrix(order,matrix)
    IMPLICIT NONE
    INTEGER i,j,order
    REAL*8,DIMENSION(order,order)::matrix
    DO i=1,order
       WRITE(*,'(10D20.9,6x)',advance='no')(matrix(i,j),j=1,order)
       WRITE(*,*)
    ENDDO
  END SUBROUTINE writematrix
END MODULE yazar

MODULE lapack
  IMPLICIT NONE
CONTAINS

  !          SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
  !*
  !*  -- LAPACK driver routine (version 3.2) --
  !*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  !*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  !*     November 2006
  !*
  !*     .. Scalar Arguments ..
  !      CHARACTER          JOBZ, UPLO
  !      INTEGER            INFO, LDA, LWORK, N
  !*     ..
  !*     .. Array Arguments ..
  !      DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * )
  !*     ..
  !*
  !*  Purpose
  !*  =======
  !*
  !*  DSYEV computes all eigenvalues and, optionally, eigenvectors of a
  !*  real symmetric matrix A.
  !*
  !*  Arguments
  !*  =========
  !*
  !*  JOBZ    (input) CHARACTER*1
  !*          = 'N':  Compute eigenvalues only;
  !*          = 'V':  Compute eigenvalues and eigenvectors.
  !*
  !*  UPLO    (input) CHARACTER*1
  !*          = 'U':  Upper triangle of A is stored;
  !*          = 'L':  Lower triangle of A is stored.
  !*
  !*  N       (input) INTEGER
  !*          The order of the matrix A.  N >= 0.
  !*
  !*  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
  !*          On entry, the symmetric matrix A.  If UPLO = 'U', the
  !*          leading N-by-N upper triangular part of A contains the
  !*          upper triangular part of the matrix A.  If UPLO = 'L',
  !*          the leading N-by-N lower triangular part of A contains
  !*          the lower triangular part of the matrix A.
  !*          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
  !*          orthonormal eigenvectors of the matrix A.
  !*          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
  !*          or the upper triangle (if UPLO='U') of A, including the
  !*          diagonal, is destroyed.
  !*
  !*  LDA     (input) INTEGER
  !*          The leading dimension of the array A.  LDA >= max(1,N).
  !*
  !*  W       (output) DOUBLE PRECISION array, dimension (N)
  !*          If INFO = 0, the eigenvalues in ascending order.
  !*
  !*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  !*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  !*
  !*  LWORK   (input) INTEGER
  !*          The length of the array WORK.  LWORK >= max(1,3*N-1).
  !*          For optimal efficiency, LWORK >= (NB+2)*N,
  !*          where NB is the blocksize for DSYTRD returned by ILAENV.
  !*
  !*          If LWORK = -1, then a workspace query is assumed; the routine
  !*          only calculates the optimal size of the WORK array, returns
  !*          this value as the first entry of the WORK array, and no error
  !*          message related to LWORK is issued by XERBLA.
  !*
  !*  INFO    (output) INTEGER
  !*          = 0:  successful exit
  !*          < 0:  if INFO = -i, the i-th argument had an illegal value
  !*          > 0:  if INFO = i, the algorithm failed to converge; i
  !*                off-diagonal elements of an intermediate tridiagonal
  !*                form did not converge to zero.
  !*
  !*  =====================================================================
  SUBROUTINE diag(order,matrix,eigenvalues)
    IMPLICIT NONE
    INTEGER :: order,info,lda
    REAL(8), DIMENSION(order,order)::matrix
    REAL(8), DIMENSION(order)::eigenvalues
    REAL(8), DIMENSION(3*order-1) :: work
    lda=3*order-1
    CALL dsyev('V','U',order,matrix,order,eigenvalues,work,lda,info)
    !    PRINT*,'info=',info
  END SUBROUTINE diag

  SUBROUTINE complex_diag(order,matrix,eigenvalues)
    IMPLICIT NONE
    INTEGER :: order,info,lda,Lwork
    COMPLEX*16, DIMENSION(order,order)::matrix
    REAL(8), DIMENSION(order)::eigenvalues
    COMPLEX*16, DIMENSION(2*order-1) :: work
    REAL(8), DIMENSION(3*order-2) :: Rwork
    lda=3*order-1
    Lwork=2*order-1
    CALL zheev('V','U',order,matrix,order,eigenvalues,work,Lwork,Rwork,info)
  END SUBROUTINE complex_diag

END MODULE lapack
